EMIT Methane Point Source Plume Complexes

Daily aggregated, global point source methane emission plume estimates from the EMIT instrument on the International Space Station (ISS)
Author

Siddharth Chaudhary, Vishal Gaur

Approach

  1. Visualize obsereved masks of the ISS and target masks of detected plumes.
  2. Identify available dates and temporal frequency of observations for the given collection using the GHGC API /stac endpoint. The collection processed in this notebook is the Earth Surface Mineral Dust Source Investigation (EMIT) methane emission plumes data product.
  3. Pass the STAC item into the raster API /stac/tilejson.json endpoint.
  4. Using folium.Map, visualize the plumes.
  5. After the visualization, perform zonal statistics for a given polygon.

About the Data

The EMIT instrument builds upon NASA’s long history of developing advanced imaging spectrometers for new science and applications. EMIT launched to the International Space Station (ISS) on July 14, 2022. The data shows high-confidence research grade methane plumes from point source emitters - updated as they are identified - in keeping with JPL Open Science and Open Data policy.

Installing the Required Libraries

Required libraries are pre-installed on the GHG Center Hub. If you need to run this notebook elsewhere, please install them with this line in a code cell:

%pip install requests, folium, rasterstats, pystac_client, pandas, matplotlib

Querying the STAC API

Please run the next cell to import the required libraries.

import requests
import folium
import folium.plugins
from folium import Map, TileLayer 
from pystac_client import Client 
import pandas as pd
import matplotlib.pyplot as plt
import branca.colormap as cm
import geopandas
from pyproj import Geod
from shapely import wkt

Description about the geojson, target mask and permian basin

# ISS data coverage 
coverage = geopandas.read_file('coverage.json')

# location where plume was detected
target_mask = geopandas.read_file('target_mask.json')

# Loading Permian Basin shape file as a region of interest
# User can pass any json or shape file here
permian_basin = geopandas.read_file('permian.zip')

# ISS area intersecting with permian basin
result_coverage = geopandas.clip(coverage, permian_basin)
# Initializing a map with a center at [43, -100] and a zoom level of 4. The use of tiles = None has been employed to remove the default basemap.

m_ = folium.Map(location=[43, -100], zoom_start=4, tiles=None)
folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png', name='ESRI World Imagery', attr='Tiles © Esri — Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',overlay='True').add_to(m_)
# Load other TileLayer from https://leaflet-extras.github.io/leaflet-providers/preview/

# From covergae.json fid, geometry has been selected
map_layer_target = folium.GeoJson(coverage[['fid', 'geometry']], name = 'ISS coverage').add_to(m_)
map_layer_target.add_to(m_)

map_layer_coverage = folium.GeoJson(target_mask, name = ' Target Mask').add_to(m_)
map_layer_coverage.add_to(m_)

# Layer control crates the toggle button for all the layers
folium.LayerControl(collapsed=False).add_to(m_)
m_
Make this Notebook Trusted to load map: File -> Trust Notebook
# Time Series Plot for observed and target mask
# Provide STAC and RASTER API endpoints
STAC_API_URL = "http://ghg.center/api/stac"
RASTER_API_URL = "https://ghg.center/api/raster"

# Please use the collection name similar to the one used in STAC collection.
# Name of the collection for methane emission plumes. 
collection_name = "emit-ch4plume-v1"
# Fetching the collection from STAC collections using appropriate endpoint.
# the 'requests' library allows a HTTP request possible
collection = requests.get(f"{STAC_API_URL}/collections/{collection_name}").json()
collection = {key: value for key, value in collection.items() if key != 'summaries'}
collection
{'id': 'emit-ch4plume-v1',
 'type': 'Collection',
 'links': [{'rel': 'items',
   'type': 'application/geo+json',
   'href': 'https://ghg.center/api/stac/collections/emit-ch4plume-v1/items'},
  {'rel': 'parent',
   'type': 'application/json',
   'href': 'https://ghg.center/api/stac/'},
  {'rel': 'root',
   'type': 'application/json',
   'href': 'https://ghg.center/api/stac/'},
  {'rel': 'self',
   'type': 'application/json',
   'href': 'https://ghg.center/api/stac/collections/emit-ch4plume-v1'}],
 'title': 'Methane Point Source Plume Complexes',
 'assets': None,
 'extent': {'spatial': {'bbox': [[-121.90662384033203,
     -39.21891784667969,
     151.0906524658203,
     50.372535705566406]]},
  'temporal': {'interval': [['2022-08-10T06:49:57+00:00',
     '2023-10-08T16:11:15+00:00']]}},
 'license': 'CC0-1.0',
 'keywords': None,
 'providers': None,
 'description': 'Methane plume complexes from point source emitters',
 'item_assets': {'ch4-plume-emissions': {'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
   'roles': ['data', 'layer'],
   'title': 'Methane Plume Complex',
   'description': 'Methane plume complexes from point source emitters.'}},
 'stac_version': '1.0.0',
 'stac_extensions': None,
 'dashboard:is_periodic': False,
 'dashboard:time_density': 'day'}

Examining the contents of our collection under the temporal variable, we note that data is available from August 2022 to May 2023. By looking at the dashboard: time density, we can see that observations are conducted daily and non-periodically (i.e., there are plumes emissions for multiple places on the same dates).

# Create a function that would search for the number of items in above data collection in the STAC API
def get_item_count(collection_id):
    count = 0
    items_url = f"{STAC_API_URL}/collections/{collection_id}/items"

    while True:
        response = requests.get(items_url)

        if not response.ok:
            print("error getting items")
            exit()

        stac = response.json()
        count += int(stac["context"].get("returned", 0))
        next = [link for link in stac["links"] if link["rel"] == "next"]

        if not next:
            break
        items_url = next[0]["href"]

    return count
# Apply the above function and check the total number of items available within the collection
number_of_items = get_item_count(collection_name)
plume_complexes = requests.get(f"{STAC_API_URL}/collections/{collection_name}/items?limit={number_of_items}").json()["features"]
parse_plume_complexes = plume_complexes
print(f"Found {len(plume_complexes)} items")
Found 752 items
# Examine the first item in the collection
plume_complex_ids = list(map(lambda d: d.get('id', f"id not found in dictionary"), plume_complexes))
plume_complex_ids[:10]
['EMIT_L2B_CH4PLM_001_20231008T161115_001520',
 'EMIT_L2B_CH4PLM_001_20231006T100206_001584',
 'EMIT_L2B_CH4PLM_001_20231006T100206_001583',
 'EMIT_L2B_CH4PLM_001_20231006T100206_001581',
 'EMIT_L2B_CH4PLM_001_20231006T100206_001580',
 'EMIT_L2B_CH4PLM_001_20231006T100206_001579',
 'EMIT_L2B_CH4PLM_001_20231006T082735_001586',
 'EMIT_L2B_CH4PLM_001_20231006T065557_001589',
 'EMIT_L2B_CH4PLM_001_20231004T174744_001453',
 'EMIT_L2B_CH4PLM_001_20231004T174744_001452']

Exploring Methane Emission Plumes (CH₄) using the Raster API

In this notebook, we will explore global methane emission plumes from point sources. We will visualize the outputs on a map using folium.

# To access the year value from each item more easily, this will let us query more explicity by year and month (e.g., 2020-02)
plume_complexes = {plume_complex["id"]: plume_complex for plume_complex in plume_complexes} 
asset_name = "ch4-plume-emissions"

Below, we are entering the minimum and maximum values to provide our upper and lower bounds in rescale_values.

# Fetching the min and max values for a specific item
rescale_values = {"max":plume_complexes[list(plume_complexes.keys())[0]]["assets"][asset_name]["raster:bands"][0]["histogram"]["max"], "min":plume_complexes[list(plume_complexes.keys())[0]]["assets"][asset_name]["raster:bands"][0]["histogram"]["min"]}

Now we will pass the item id, collection name, and rescaling_factor to the Raster API endpoint. We will do this for only one item so that we can visualize the event.

color_map = "magma" # please refer to matplotlib library if you'd prefer choosing a different color ramp.
# For more information on Colormaps in Matplotlib, please visit https://matplotlib.org/stable/users/explain/colors/colormaps.html

# Select the item ID which you want to visualize. Item ID is in the format yyyymmdd followed by the timestamp. This ID can be extracted from the COG name as well.
plume_complex_id = "EMIT_L2B_CH4PLM_001_20230418T200118_000829"
methane_plume_tile = requests.get(
    f"{RASTER_API_URL}/stac/tilejson.json?collection={plume_complexes[plume_complex_id]['collection']}&item={plume_complexes[plume_complex_id]['id']}"
    f"&assets={asset_name}"
    f"&color_formula=gamma+r+1.05&colormap_name={color_map}"
    f"&rescale={rescale_values['min']},{rescale_values['max']}", 
).json()
methane_plume_tile
{'tilejson': '2.2.0',
 'version': '1.0.0',
 'scheme': 'xyz',
 'tiles': ['https://ghg.center/api/raster/stac/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?collection=emit-ch4plume-v1&item=EMIT_L2B_CH4PLM_001_20230418T200118_000829&assets=ch4-plume-emissions&color_formula=gamma+r+1.05&colormap_name=magma&rescale=-638.1588745117188%2C2034.2767333984375'],
 'minzoom': 0,
 'maxzoom': 24,
 'bounds': [-104.76285251117253,
  39.85322425220504,
  -104.74658553556483,
  39.86515336765068],
 'center': [-104.75471902336868, 39.85918880992786, 0]}

Visualizing CH₄ Emission Plume

colormap = "magma" # please refer to matplotlib library if you'd prefer choosing a different color ramp.
# For more information on Colormaps in Matplotlib, please visit https://matplotlib.org/stable/users/explain/colors/colormaps.html

#Defining the breaks in the colormap 
color_map = cm.LinearColormap(colors = ['#310597', '#4C02A1', '#6600A7', '#7E03A8', '#9511A1', '#AA2395', '#BC3587', '#CC4778', '#DA5A6A', '#E66C5C', '#F0804E', '#F89540','#FDAC33', '#FDC527', '#F8DF25'], vmin = 0, vmax = 1500 )
color_map.caption = 'ppm-m'

# Select the item ID which you want to visualize. Item ID is in the format yyyymmdd followed by the timestamp. This ID can be extracted from the COG name as well.
plume_complex_id = "EMIT_L2B_CH4PLM_001_20230418T200118_000829"

methane_plume_tile = requests.get(
    f"{RASTER_API_URL}/stac/tilejson.json?collection={plume_complexes[plume_complex_id]['collection']}&item={plume_complexes[plume_complex_id]['id']}"
    f"&assets={asset_name}"
    f"&color_formula=gamma+r+1.05&colormap_name={colormap}"
    f"&rescale={rescale_values['min']},{rescale_values['max']}", 
).json()
methane_plume_tile

# Set initial zoom and center of map for plume Layer
map_ = folium.Map(location=(methane_plume_tile["center"][1], methane_plume_tile["center"][0]), zoom_start=13, tiles=None, tooltip = 'test tool tip')
folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png', name='ESRI World Imagery', attr='Tiles © Esri — Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',overlay='True').add_to(map_)


map_layer = TileLayer(
    tiles=methane_plume_tile["tiles"][0],
    name='Plume Complex',
    overlay='True',
    attr="GHG",
    opacity=1,
)
map_layer.add_to(map_)
folium.LayerControl(collapsed=False, position='bottomleft').add_to(map_)
map_.add_child(color_map)
# visualising the map
map_
Make this Notebook Trusted to load map: File -> Trust Notebook
# Creating a dataframe for each plume in the STAC and its corresponding area 

plume_df = pd.DataFrame()
for plume in parse_plume_complexes:
    temp_polygon = geopandas.GeoDataFrame.from_features([plume])['geometry'].values[0]
    geod = Geod(ellps='WGS84')
    ply = wkt.loads(str(temp_polygon))

    temp_dict = {'id':plume['id'], 'geometry':temp_polygon, 'area(km2)':(abs(geod.geometry_area_perimeter(ply)[0])/1e+6)}
    plume_df = plume_df._append(temp_dict, ignore_index = True)
plume_df
id geometry area
0 EMIT_L2B_CH4PLM_001_20231008T161115_001520 POLYGON ((-103.94950373078798 31.8037824889992... 0.605059
1 EMIT_L2B_CH4PLM_001_20231006T100206_001584 POLYGON ((3.2823994751112684 33.47982320911174... 6.503630
2 EMIT_L2B_CH4PLM_001_20231006T100206_001583 POLYGON ((3.4271755580197185 32.97608919779358... 10.971985
3 EMIT_L2B_CH4PLM_001_20231006T100206_001581 POLYGON ((3.4889900653289443 33.00320082380639... 5.484339
4 EMIT_L2B_CH4PLM_001_20231006T100206_001580 POLYGON ((3.5101371336189424 33.04224156526485... 9.456115
... ... ... ...
747 EMIT_L2B_CH4PLM_001_20220811T042630_000490 POLYGON ((56.39299036918186 22.31688231459388,... 12.543665
748 EMIT_L2B_CH4PLM_001_20220810T065132_000496 POLYGON ((40.48659938746134 35.26702159587666,... 11.778459
749 EMIT_L2B_CH4PLM_001_20220810T065132_000487 POLYGON ((40.442136320800316 35.29250652432872... 6.625069
750 EMIT_L2B_CH4PLM_001_20220810T065021_000486 POLYGON ((36.22139838312476 31.89813094752386,... 2.590405
751 EMIT_L2B_CH4PLM_001_20220810T064957_000485 POLYGON ((35.192241059678175 31.11189379315212... 7.524315

752 rows × 3 columns

Visualizing CH4 Emission Plumes in the Permian Basin

# Converting plume_df into a geo dataframe for further geospatial analysis
geo_temp_df = geopandas.GeoDataFrame(plume_df, geometry=plume_df.geometry, crs = 'EPSG:4326')
# Selecting all the plume tagets which are under permian basin
geo_temp_df = geopandas.clip(geo_temp_df, permian_basin)
colormap = "magma" # please refer to matplotlib library if you'd prefer choosing a different color ramp.
# For more information on Colormaps in Matplotlib, please visit https://matplotlib.org/stable/users/explain/colors/colormaps.html

color_map = cm.LinearColormap(colors = ['#310597', '#4C02A1', '#6600A7', '#7E03A8', '#9511A1', '#AA2395', '#BC3587', '#CC4778', '#DA5A6A', '#E66C5C', '#F0804E', '#F89540','#FDAC33', '#FDC527', '#F8DF25'], vmin = 0, vmax = 1500 )
color_map.caption = 'ppm-m'
# Select the item ID which you want to visualize. Item ID is in the format yyyymmdd followed by the timestamp. This ID can be extracted from the COG name as well.

map_ = folium.Map(location=[32,-102], zoom_start=7.4, tiles=None)
folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png', name='ESRI World Imagery', attr='Tiles © Esri — Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',overlay='True').add_to(map_)

for tile_id in geo_temp_df['id']: 
    methane_plume_tile = requests.get(
        f"{RASTER_API_URL}/stac/tilejson.json?collection={plume_complexes[tile_id]['collection']}&item={plume_complexes[tile_id]['id']}"
        f"&assets={asset_name}"
        f"&color_formula=gamma+r+1.05&colormap_name={colormap}"
        f"&rescale={rescale_values['min']},{rescale_values['max']}", 
    ).json()
    methane_plume_tile

    map_layer = TileLayer(
        tiles=methane_plume_tile["tiles"][0],
        name=tile_id,
        overlay='True',
        attr="GHG",
        opacity=1,
    )
    map_layer.add_to(map_)

map_layer_permian = folium.GeoJson(permian_basin, name= 'Permian Shape').add_to(m_)
map_layer_permian.add_to(map_)
folium.LayerControl(collapsed=True, position='bottomleft').add_to(map_)
map_.add_child(color_map)
# visualising the map
map_

# Zoom in to visualize different plumes
Make this Notebook Trusted to load map: File -> Trust Notebook

Summary

In this notebook we have successfully explored, analyzed, and visualized the STAC collection for EMIT methane emission plumes.